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ABSTRACT 



In most references, for the great convenience in calculations, the exterior gravity potential of 
the spheroidal earth is presented as a series of spherical harmonic terms. In an early part of this 
paper that series is derived in full detail as a classical solution of Laplace’s equation in spherical 
coordinates. In a later part of the report the series is truncated at the l/r 3 terms as is often done 
in the literature. It is shown that this truncation gives the potential for a rotating triaxial earth, 
having small oblateness and equatorial eccentricity, with an equipotential surface. The coefficients 
of the truncated series are interpreted as functions of the geometric parameters of the matching 
triaxial earth model and its rotation. The harmonic coefficients, developed according to the theory 
presented in this paper and measured to date from gravity surveys and satellite data, show that the 
earth has: 

1. an oblateness of 1/298.2 

2. an equatorial eccentricity of 1.14 x 1(T S , which corresponds to a difference of major and 
minor equatorial diameters for the earth of about 475 feet. 

In addition to the harmonics contributing principally to the triaxial earth, the higher order 
harmonics are interpreted in the latter part of the paper as reflecting additional small longitude 
and latitude-dependent surface deviations from this triaxial earth. The effect on the earth surface 
model (the geoid) of terms through J 44 is calculated, but the longitude-dependent deviations them- 
selves are left in general terms as functions of the harmonic coefficients. These latter coef- 
ficients must be regarded as largely unknown except for the second tesseral harmonic J 22 which 
(together with J 20 and the centrifugal potential) fixes the surface shape of the triaxial geoid. Longi- 
tude independent (zonal) geoid deviations through j 40 are calculated from a zonal potential due to 
Kozai (1962), which is probably the most accurate to date. 

While the details of the topography of the geoid remain to be discovered by geodetic research, 
certain large-scale features appear to be discernable already. These are discussed in a later 
section of this report. Except perhaps for J 22 determined from 24 hour satellite data, it appears 
necessary at this date to consider all longitude gravity effects as a set through at least the fourth 
order in the long term applications. 

Only the second order longitude gravity effect probably needs to be considered in the appli- 
cations concerning the 24-hour satellite, thou gh the effects of third- and fourth-order earth gravity 
may be discernable over long periods of time. 

In the last section of this report the gravity and geometric effects of the small nutations of the 
north pole are discussed. The gravity effects are shown to be entirely negligible for most applica- 
tions. The geometric effects might well be discernable by close long-term study of the geoid 
coordinate errors of stations used as references in the observation of earth satellites. 
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THE GRAVITY POTENTIAL 
OF THE EARTH THROUGH FOURTH ORDER 


by 

C. A. Wagner 

Goddard Space Flight Center 


I. INTRODUCTION 

The basic purpose of this report is to gather in one document information on the earth’s ex- 
ternal gravity field necessary for the accurate determination of spacecraft trajectories in its 
vicinity. Much of the material on the mathematical representation of the field is not new. In most 
references, for convenience of calculation and because the earth itself is so nearly spherically 
symmetric, the general gravity field has been developed into an infinite series of spherical har- 
monic components. The first two terms of this series have been known to some accuracy since 
Isaac Newton’s time. In recent years, with the results of close tracking of artificial satellites 
available, many more terms of this series have become known. It seemed strange, with so much 
history behind it, that in 1962 when the author began his investigations of the gravity effects on a 
24-hour satellite, no single reference could answer adequately the basic question: Where does the 
familiar harmonic series representation of the earth’s gravity field come from? That one question, 
and the attempt to answer it, gave rise to this document. But in the course of answering that basic 
question other questions arose on the physical and geometric interpretation of the individual gravity 
components for which, likewise, no adequate answers could be found in the literature. The inter- 
pretation of satellite observations is always with respect to specific stations fixed on the earth’s 
surface. The form of that surface, called the geoid or mean sea level surface, is a reflection of 
the mass distribution of the earth beneath it just as the gravity series is. Thus the problems of 
accurate orbit determination and accurate ground station locations are interrelated. Their 
common point of departure is the geometric (geoid) interpretation of the terms in the gravity 
field. These are just some of the additional questions about the field which are important in 
the applications to geodesy and satellite orbit determination and for which answers are given 
in this report. 

For example, in the last section of this report the gravity and geometric effects of the small 
nutations of the north pole are discussed and calculated, so far as is known, for the first time in 
the literature. 
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II. THE GRAVITY POTENTIAL AND THE GRAVITY FIELD OF THE EARTH 


In this study we are interested fundamentally in determining the exterior gravity force field of 
the earth in a form most efficient for the calculation of the trajectories of small bodies in the 
vicinity of the earth. In Figure 1 the earth is represented by an extended domain, D e , in the space 

Xj, X 2 , X 3 . The gravitational force exerted by 
the small mass dm in the extended earth domain 
D e on a particle of unit mass at P is, from 
Newton’s law of gravitation, 



K), 


- Gdm ^ 

P 


(i) 


= - Gdm p 
on a unit p 2 

mass at P 


where G is the universal gravity constant 
(= 6.67 x 10~ 8 dynes-cm 2 /gm 2 ). The total 
gravity force of D e on the unit mass at P is the 
integral of Equation 1 over D e , 


Figure 1— The gravity force field of the earth. 

The unit vector p can be decomposed as follows: 


(fJ, 


I 


-Gdm 


(2) 


P _ 


x i (x, -fj) +s 2 (x 2 -£ 2 )+S 3 (xj-^ 3 ) 


1 1/2 


{( X l _ ^l) 2 + ( X 2 “^2 ) 2 + ( X 3 _ ^ 3 ) 2 } 

Thus Equation 2 can be rewritten as the vector equation 

{M X i“^i) +X 2 ( X 2"£ 2 ) +X 3 ( x 3"^3)} 


F. = - G 


I 


dm 


. {( X l"^l) 2+ ( X 2'^ 2 ) 2+ ( X 3”^3) 2 } 3/2 


(3) 


Equation 3 implies that the three scalar equations for the (F e ) x , (F e ) x and (f c ) x coordinate com- 
ponents of total attraction, (F e ) p = (F e ) Xj x t + (F e ) x 2 + (F e ) Xl x 3 , are 

( x i 


( F ) = - G f , 

‘ k {( X l-^l) 2 + ( X 2-^2) 2+ ( X 3-^3) 2 } 

_ (• (x, “ ) dm 

" ' 0 J . — 7 — 


i = 1 , 2 , 3 ■ 


(4) 
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For an arbitrary mass distribution in D e it would be extremely cumbersome to compute the gravity 
field through the three triple integrals of Equation 4. For most applications it is far more desir- 
able and instructive to compute the field with a single intermediate scalar function V, called the 
potential, whose easily computed gradient (w = x t dv/dx 1 + x 2 dV/'dx 2 + x 3 d\/dx 3 ) gives the force 
field represented by Equation 4. The gravity potential of D e at P, (V e ) p is defined by the scalar 
triple integral 


( V e). 


f Grim 

I — = G 


f rim 


( 5 ) 


Note this implies that each small mass contributes d(v e ) p = Ck Wp to the total potential. For 
example, the scalar potential of a point mass, m, on a unit mass at p from m in space is 


The gradient is a variation of the field with respect to the coordinates of the field point only 
and is therefore independent of the limits of integration for the domain D e . Thus the V operation 
can be taken on the integrand and the gradient of Equation 5 with respect to the space coordinates 
of the arbitrary field point P is 


But 


V ( V e)r = V f ^ = G J dmV (^) 

J D Jd 


and 




('j '«•)’] [ 2 ( x i -*.)]} 


P 


i = 1, 2, 3 . 


Therefore, from Equation 8 and the definition of the gradient Equation 7 becomes: 


' 7 ( V e )p ~ G p 3 { X l( X l ^l) +X 2 C X 2 ^2 ) X 3 ( X 3 ^3)} 


(V 


( 8 ) 


3 



\ 


or 


*( V e)p „ ^ ( V e) P . a(V e ) 

+ X 


^Xj T x 2 (?x 2 + X 3 (9x 3 


hi 

D 


i, l-G 


( x i 


+ 


" G I 

•'D 


(x 2 -^ 2 ) dm' 


+ X, 


- G 


1 


(x 3 -^3) dm" 


( 9 ) 


Equation 9 is equivalent to three scalar equations: 


d ( V e) 

(9x . 


f ( x i-^' 

J D p 3 


p _ 1 ( x i _ ^i) 

- G 


= 1, 2, 3 . 


(10) 


From Equations 10 and 4 the components of the gravity force and gravity potential gradient vectors 
are identical: 


<H v .)p 

dx. 

i 



i = 1, 2, 3 • 


(id 


Thus the problem of finding the vector gravity field of the earth has been broken down to the simpler 
problem of finding the single scalar gravity potential of the earth (Equation 5). The three gradient 
derivatives of this scalar function will be the coordinate components of the vector gravity force 
field (Equation 11). 

III. THE EXTERNAL POTENTIAL OF THE EARTH 
AS AN INFINITE SERIES OF MASS INTEGRALS 

In determining the gravity potential of the earth at P from Equation 5, there are three methods 
of proceeding: 

1. If the precise nature or form of the mass domain of the earth is known or can be assumed, 
the elemental potential Gdm/p can be expressed in accordance with this form to allow an exact 
integration of Equation 5. In such a manner Newton first proved that the external gravity potential 
of a spherical mass whose density is only a function of distance from the center of mass (c.m.) is 
given simply as 


( V) p - CM/r , (12) 

where Mis the total mass of the extended sphere and r is the distance from the c.m. of the 
sphere to the field point P. In other words, Equation 12 states that an extended earth of 
radially symmetric mass distribution can be treated as a mass point with a simple central 
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foroe field given by 


( F e)l 



(13) 


Since the real mass distribution of the earth can be closely approximated by such a spherically 
symmetric mass domain, the exact form of the earth's potential is Equation 12 with small, but 
important departures which depend, in general, on all three coordinates of the field point inde- 
pendently and not just the distance from the c.m. of the earth. 

2. In the general case of a nearly spherical mass domain the integrand of Equation 5 can be 
expanded in an infinite harmonic series of terms in lA n (n = 1, 2, 3, * * •). Integration of the leading 
term (l/r) will give the dominant spherically symmetric part of the potential. The remaining 
series, integrated term-by-term with respect to the domain D e , will give mass integrals of small 
order related to the nonspherical part of the mass distribution. 

3. In the general case of an arbitrary domain the indicated mass integration over the domain 
in Equation 5 for V e can be replaced by a differential equation which V e must satisfy independently 
of the nature of the domain. The solution of the differential equation in spherical coordinates as 
an infinite series of spherical harmonics will yield the same series solution as in method 2. But 
in this case, the coefficients of the series will be arbitrary. Spherical coordinates are chosen to 
give a series whose first few terms approximate well the nearly spherically symmetric earth. The 
form of this series is also convenient for calculating the perturbations of nearly circular satellite 
orbits, as well as displaying the gravity field on the earth’s nearly spherical surface. The coef- 
ficients themselves can be related to the integrated mass properties of the earth by matching 
term-by-term with the solution from method 2 (see part VI). A large number of coefficients of 
this series have been determined to date by many kinds of geodetic measurements which essentially 
map and sample the gravity field over the earth’s surface (gravimeter surveys and astro -geodetic 
survey mappings) and in space (satellite perturbation studies). 

Proceeding with the direct integration of Equation 5 (see Figure 1), we note from the law of 
cosines that 


g2 + r 2 - 2 £t cos \p , 


or 

p - r {l + (f/r) 2 -2(£/r) cos^} 1/2 
Equation 14 in Equation 5 gives (noting r is independent of D e ) 


(14) 


( 


V e)p = 7 I dm [l + (^/r) 2 -2(^/r) cosi/-]' 1/2 




1/2 


1 - 


Mlf 


- 1/2 


(15) 
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It is evident from Equation 15 that for field points far from the earth (assuming the center of co- 
ordinates is somewhere inside the earth), (V e ) p far from the earth - GM e /r. This implies the earth acts 
as an effective point mass (independent of the exact nature of its mass distribution) at a point suf- 
ficiently far from its vicinity. 

For exterior points P for which (£/r) < 1, since 2(£/r) cos<///l + (£/r) 2 < 1, the integrand of 
Equation 15 can be expanded into two convergent binomial series, yielding 


(v.)p - ^ 


dm y - — 2 — + 


•j"^ 1 * (t) c<>> * [*' (4) *•••]* -T^- CO*’# [l- j(f) *•"] 




Collecting terms of Equation 15a in powers of (£/r) gives 


( V .)p r dm <|l+ (y)cosi//+ [2-3sin 2 ^] + o(-|) + • • • ^ 


If the origin of coordinates is the c.m. of the earth, 


I dm[^ cos 0] 

J D 


0 , 


- (15a) 


(16) 


(17) 


from the definition of the c.m. of a body. This property is independent of the kind of body-fixed 
coordinate system chosen. Thus, for orthogonal axes Xj 1 , x^, x* with origin at the c.m. and x* 
oriented along r to the field point P, = fcos^, and f dm(£ cos \p) = j g* dm = 0 . Also note 
from Figure 1 that 

I £ 2 dm = I„ , (18) 

J D 

e 


the total moment of inertia of the mass of the earth about its center of mass (c.m.). This quantity 
is independent of the field point P. In addition note that 


dm ( £ s i n i// ) 2 - I f , ( 19 ) 

-V 

the moment of inertia of the earth about r. This quantity is a function of the direction of P from 
the c.m. of the earth, but not its distance. Integration of Equation 16 term-by-term with the 
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results of Equations 17, 18 and 19 gives the first two mass integrals of the external gravity poten- 
tial of the earth as 

(M 21 n - 31 / \4 1 

\ e 0 r . l 


( y .)p = * •}• < 2 °) 

l r can be shown to be a simple function of the direction of P from the c.m. and the three principal 
moments of inertia of the earth. 

For example, let x t , x 2 , x 3 be principal c.m. fixed axes for the earth and A, B, c be the cor- 
responding principal moments of inertia for the earth defined as: 


[ dm (^2 2+ ^3 2 ) 

e 

[ dm ( 5 , 2 + <f 3 2 ) 
•'D 


But from the law of cosines 


f dm(^ , 2 + € 2 ) 

Jd 


[ = I £ 2 sin 2 \pdm -*• I £ 2 (l - cos 2 \jj) dm - I dm c 2 1 

J D Jb 


Ujri 2 - p 2 ) 2 


r f [f, 2+ f 2 2+ ^3 2+X l 2+X 2 2+X 3 2 -( X l _ ^l) 2 ‘( X 2 _ -2) 2_ ( X 3*^3) 2 ]‘ 

' r Jd 4(x, 2 +x 2 2 +X 3 2 ) 


J dm ^ 

J D L 


2 + S , 2 + f , 2 - 


(x, g, +x 2 ^ 2 +x 3 ^ 3 ) : 
X l 2 + x 2 2 + X 3 2 


X l 2 (^2 2 + ^3 2 ) +X 2 2 Ul 2 + ^3 2 ) +X 3 2 (^1 2 + ^2 2 ) ~ 2 ( X 1 x 2^l ^2 +X 1 X 3 ^1 js + X 2 X 3 *2 fjjl 


But, for principal axes of inertia Xj , x 2 , x 3 , 


and therefore 


j fj^jdm - fi£ 3 dm - i 2 £ 3 

Jo* J De J D e 

I r = [A(x 1 ) 2 +B(x 2 f +C(x 3 )*]/r 


dm - 0 


M 1 + Brn 2 + Cn 2 
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* 

where m and n are the direction cosines of P from the c.m. of the earth with respect to earth- 
fixed principal axes of inertia. From Equation 18 we have: 


2l r 


4 = 4 [£ (£,’+#/) (fi 2+ ^ a ) + l(^ 2+ 0]*' 

Jd 


A + B + C 


(20c) 


For a spherically symmetric earth mass distribution 


A = B = C , 

and from Equations 20c and 20b, 2I 0 = 3A and 

31, = 3 [a(xj) 2 + B(x 2 ) 2 +C(x 3 ) J ] /r 2 = 3A [ (x,) 2 + (x 2 )’ + (x,)’ ]/r 2 = 3A (20d) 

(for a spherically symmetric earth). Thus from the above results in Equation 20 it is seen that the 
1/r 3 term of the gravity potential is zero for a spherically symmetric mass distribution. Thus, for 
a triaxial spheroidal earth, not differing greatly from a spherically symmetric mass domain, the 
second term of the potential function, Equation 20, will be close to zero even for field points in the 
near vicinity of the earth. The same argument can be made for the higher order terms of Equation 
20. These can be shown to measure earth mass deviations with planes of symmetry not necessarily 
orthogonal. 

A summary of the development of the gravity potential of the earth thus far shows: 

1. The first term (the 1/r term) of the infinite harmonic series of mass integrals of the 
gravity potential of the earth gives the effect of the earth 1 s mass distribution concentrated at a 
point. 

2. The second term (the l/r 3 term) gives the gravity effect due to the deviation of the earth’s 
ellipsoid of inertia from a sphere. The 1/r 2 term in the potential is absent if the c.m. is the center 
of the coordinate reference. 


IV. THE EXTERNAL GRAVITY POTENTIAL OF THE EARTH 
IN TERMS OF LAPLACE’S EQUATION 

The Gravity Potential as a Solution to Laplace’s Equation 

The operator V 2 ( ) of a scalar function ( ) is defined as 
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Taking V 2 (vj p in Equation 5 and noting that this operation is independent of the manner in which 
the mass domain D e is integrated gives 


V2 (V e )e 




dm V 2 (1/ p) 


(21) 


But 


and 


72 (1/p) 


<? 2 (1/p) + d 2 (1/p) + <? 2 (1/p) 

d x 2 d x 2 5 x 3 2 


a 2 (1 /P) 
a x 2 

i 


s (^) 

(9x i <9x. 


dx i 

<2x. 


Thus 


- p 3 + 3(x. -£.) 2 p 5 ’ * 1* 2. 3 . 

V 2 (l/p) = -3/T 3 + 3p' s |( X j ~£,) 2 + ( x 2 -£ 2 ) 2 +( x 3 -£ 3 ) 2 } 


= -3p' 3 + 3p~ s P 2 = 0 . 

Thus since the integrand in Equation 21 is zero at every mass point (<fj , Z 2 , <f 3 ) in the domain D e , 
the integral is also zero. Therefore, the exterior gravity potential (where p > 0 everywhere, and 
assuming no singularities in the potential at P) must satisfy Laplace's Equation, 


? 2 (V.)p = 0 . (22) 

(See Reference 1, Section n for a more complete discussion of the derivation of Equation 22 and the 
mass integral potential series. Equation 20.) 
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The External Gravity Potential as an Infinite Series of Spherical Harmonics 


From the suggestion of part HI, since the earth is nearly spherically symmetric in form, 
solutions of Laplace 1 s equation in spherical coordinates will be developed; and the number of terms 

necessary to obtain an accurate gravity field 
x, for most purposes will be a minimum. (El- 

lipsoidal coordinates would require fewer terms 
for accuracy but they would be computationally 
more cumbersome.) 

Consider body "D e " with the origin of 
coordinates x % , x 2 , x 3 fixed inside "D e ", and 
the exterior field point p specified by spherical 
coordinates <£, L, r (Figure 2). 

In general orthogonal coordinates q t , q 2 , q 3 



V 2 V 


1 

Ll 

Q 2 ^3 dV 

<?> Q 2 Q 3 ' 

1 

Qi 


dq 7 


Qj Q 3 6 >v 

Q 2 dq 2 




Q 1 Q 2 dV 
Q, dq. 


(23) 


where ds £ r q a dq i , a small line element in the direction of increasing q i coordinate at P, all other 
coordinates being constant (Reference 2, pp. 173 and 175). 

From the spherical coordinates of Figure 2 the following terms are identified: 


ds 4, ~ 

r cos Ld<^, 

implying - r cos L 

ds L = 

r dL , 

implying Q L n r 

ds ~ 

r 

1 dr , 

implying Q r = 1 . 


With these identifications in Equation 23, Equation 22 becomes 


l 

r 2 cos L 


/ 9 

r 1 avi 

x d 

r T c> v i 

d 

r , <?vi\ 


_cos L d$} 

+ 1C 

[ C °S L dh 

+ dr 

L r cosL^JJ 


Expansion of this equation will give 


1 

i a 2 v 

d 

r . iYi 

r 2 cos L 

l cos L d( f,2 

dL 

L cos L dc J 


d 

cosL^r 
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or 

4 


1 d I 2 1 i I , 1X\ + __L_ ill - „ 

r 2 ' 9 r ) r 2 cos L ^ ' ° S 9h J r 2 cos 2 L d4> 2 ' 


( 24 ) 


To separate variables, assume 


V = R(r)H(L )<>(<£) . 


(25) 


Substituting Equation 25 into Equation 24, taking the indicated partial derivatives, and rearranging 
terms will transform Equation 24 to 


m> a 

t 2 


a /. a r\ gj> a t , sh\ rh a 2 <i> 

9 r V 9 t I + r 2 cos L 9L \ cos L aL / + r 2 cos 2 L d< ^2 

Multiply Equation 26 by r 2 cos 2 L to get 

, a / , aR\ _ a ( , aH\ a 2 <t> 

»»cos 2 L^r (r 2 ^r) + R4» cos L ^ (cosL^j+RH — 


0 . 


= 0 


(26) 


(27) 


Next divide Equation 27 through by Rdh to isolate <J>(-i>) . Equation 27 now becomes 


cos 2 l a / , aR \ cos l a / , aH\ 

r aT \ r aF/ + ~ i h“ aE V cos L aL ) 


l a 2 <e 

"• a<t> 2 


+ 


(28) 


Small m is a constant (to be determined later) since the far left-hand side of Equation 28 is a 
function of r and L only, while the right-hand side, (- l/<i> d 2 ) , is a function of </> only. 


From Equation 28 


d 2 Q 

d<f > 2 


m 2 $ 


0 , 


which implies 


<t>(4>) = K 0 cos m(<£ ~<P 0 ) , (29) 

where K 0 and <p 0 are constants of integration. Since V must be periodic at least every 2 77 radians, 
and cosm(<£-c£ 0 ) = cos , it is necessary and sufficient that m be restricted to all positive 

integers including zero. 
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Dividing Equation 28 by cos 2 L and rearranging terms allows one to isolate H and R. Equation 28 
now becomes 


1 d( 8H\ m 2 _ 1 d ( , dR\ _ 

hToTl dh \ cosL It) ~ - "If 57 V 17) - ~ q ■ (30) 

where q is a constant (to be determined later), since the far left-hand side of Equation 30 is a 
function of L only and the near right is a function of r only. 

From Equation 30 


_LA 

R Sr 



q . 


and also 


i .a / _ *h\ 

H cos L dL \ cos L dL I 


cos 2 L 


q . 


(31) 


(32) 


The independent variable in Equation 32 is transformed from L to x by setting x - sinL. Now 

d( ) _ d( ) dx d( ) . 

dL dx dL cos ^ dx ’ 

and 

cos L n - x 2 ) 1/2 , cos 2 L = 1 - x 2 , etc. (33) 

In Equations 31 and 32 let q n n(n + 1), where n is also a constant which will be determined later. 
With these changes Equation 32 becomes 


JL A 

H dx 


(i 


x ) dx 



+ n(n + 1 ) 


0 . 


(34) 


Note that the partial differentiation in Equation 32 is equivalent to total differentiation since only 
one independent variable (L) is involved. 


Expanding Equation 34 and multiplying through by H gives 


r d 2 H ^ dH 

[l x ^ dx 2 2x dx + 


n(n + 1 ) 


1 -x 2 J 


(35) 
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Equation 35 is recognized as the associate Legendre equation . It has applicable solutions 


H 


nm 


DnJl -* 2 )" 172 


d"P n (X) 
dx m 


(Reference 3, 

pp. 311-313) 


( 36 ) 


The P n (x) are the well known Legendre polynomials (Reference 3, pp. 305, 306, 312 and 313), and 
m and n are non-negative integers with m < n . H nm refers to the harmonic of order n and power m . 

Other solutions of 35 are possible when n is or is not an integer if | x | / 1. But |x| = | sin L| < l 
for the exterior space suiruunding the earth, and the other solutions are thus nut admissible. They 
give infinities for |x| =1 (Reference 4, pp. 192-194 and 199). 

The polynomials P n (x) take the form 

p (X ) = _1ML/ X n . n(n-l)x n ' 2 | n(n- l)(n-2) (n-3)x n ~ 4 _ . 1 . ( 37 ) 

2 " (n! ) 2 \ ( 2 ')(l!)( 2 i-l) ( 2 a ) (2 , ) (2n - 1 ) ( 2 n - 3 ) J 

It may be noted from Equation 37 that if n is not an integer, the polynomial series becomes infinite 
and complex and thus is not applicable to the solution for V . A full treatment of the derivation of 
the Associated Legendre Functions H can be found in Reference 13. The H are there shown to 

n m nm 

be a complete, orthogonal set of functions giving the only finite solutions to Equation 35 for |x| - 1. 
In particular, the first six polynomials are: 


Po 

(X) 

- P 0 (sinL) 

=• 1 


Pi 

<*) 

- (sinL) 

= sinL 


p 2 

(X) 

- P 2 (sin L) 

= ^ (3 sin 2 L - 1 ) 


P 3 

(X) 

= P 3 (sin L) 

- ~2 (5 sin 3 L _ 3 sin L ) 


?4 

(X) 

~ P 4 (sin L) 

= (1/8) (35 sin 4 L- 30 sin 2 

L + 3) 

P 5 

(X) 

= P 5 (sin L) 

- (1/8) (63 sin 5 L + 70 sin 3 

L + 15 sin l) 

P 6 

(X) 

= P 6 (sin L) 

= (1/16) (231 sin 6 L“ 315 s 

in 4 L + 105 sin 2 L “ 5 ) 

To determine the 

variation in r 

, in Equation 31 we set q = 

n(n + 1 ) and expand to get 


<? 2 R 


<3R 


^7 r + 2r Jr ' n ( n + !) - 0 • 

Equation 39 is a "Cauchy" equation (Reference 4, p. 199) whose general solution is 


(39) 


R(r) = A, r" + A 2 r~< n+1 > . 


(40) 
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If (as is conventional) we arbitrarily stipulate that V e -0 as r -<», then Aj must be set equal to zero 
to avoid the infinities from the first term in Equation 40 for n > o and to avoid a finite v e as r 
for n - o . 

When one combines Equations 29, 36 and 40 and formulates new coefficients F nm and <p nm9 one 
obtains the series solution of Laplace’s equation in sperhical coordinates (for V e ) as 


v 


LIj- 

n-0 n>~0 


, d m P n (sinL) 

r"( n+1 ) cos™ L 

d( sin L) m 


cos m 


(*-0 


(41) 


The potential field of Equation 41 is usually separated into zonal (m - o) and sectorial-tesseral 
(or longitude) terms (m f o): 


\ — 1 X — ' V” 1 d m P n (sinL) 

Ve = 2^ F "0 r ‘ <n + ,)P n ( SinL > + ••• + r ' <n+1) COS " L CO* - (42) 


The sectorial terms of Equation 42 are those of m / o for which n = m. The tesseral terms (dis- 
tinguished from the sectorial terms) are those of m / o for which n f m . The 0 nrn represent the 
location in the space ( 0 , L, r ) of the principal plane of longitudinal mass symmetry for the nm 
harmonic (m / 0 ). 


V. MATCHING THE FIRST TWO MASS INTEGRALS OF THE POTENTIAL 
WITH THEIR EQUIVALENT SPHERICAL HARMONIC REPRESENTATIONS 

The series in Equation 42 terminated at n = 2 will give 

V e = T F 00 + [ F 10 SinL + F ll C0SLC0S + 7^ [ 3F 22 COs2LCOS 2(^6 * 022 ) 


F 20 . 3Fj j 1 / c \4 

+ s i n2 L ~ 1 ) + — sin 2L cos (0 - 0 21 ) + order [~j + * * * * (43) 

Consider the potential of a finite continuous distribution of matter at a point far enough from 
the distribution that the radius from the origin of coordinates to the field point is greater than the 
radius to any part of the distribution. If the origin of coordinates is taken at the center of mass of 
the distribution, it was shown (in Section HI) that the expansion (l/p) in Equation 15a is absolutely 
and uniformly convergent, which justifies the term-by-term integration in Equation 20 to give 
the result 



Me + ( 2 Iq - 3 I r ) 


order 



+ 
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« The discussion of the physical significance of the coefficients of Equation 43 is simplified by 
considering axes X,, X 2 , X 3 corresponding to the principal moments of inertia A, B and C of the 
mass distribution (Figure 2). 

For the field represented by the mass-integrals of Equation 20, the expansion in harmonics of 
Equation 43 is equivalent to the mass-integral harmonic series of Equation 20 if and only if 


Foo = GMe 


(44) 


b [ F 10 sinL + Fji cosLc °s (</>-<*>!,)] = 0, 


(45) 


3F, 


A 2 0 / . _ \ or 2 1 r / G ( 2 I 0 - 3 I r) 

~Y (3 sin 2 L ~ 1 ) + ~ 2 “ |_ sin X cos (0“0 2 i)J + 3F 22 cos 2 Lcos 2{<p~<P 22 ) = 2 ’ etc - (46) 


Equation 44 gives the familiar GM e /r potential term which, if it were the only term of the series 
present, would state that the whole mass of the earth, M e , acted effectively at its c.m. Equation 45 
implies that F 10 = F n = 0 since the two terms are linearly independent. The following new constants 
are defined 


( F 2 l)o 

3F 2 i 

- 2 Q cos 4> n , 

( F 22)o 

F 22 

- “G“ COS 2022 > 


( F 2 o)o 

t?\° 

II 

( F 22) 1 

F 

r 22 

- -Q- sin 20 22 . 

(46a) 

( F 2l)l 

3F 21 . 

~ 2G sin ^2i' 





To facilitate the term-by-term matching of Equation 46, transform back to rectangular principal 
coordinates x x , x 2 , x 3 , where 


Xj = r cosL cos 0 , 


x 2 - r cos L sin <p 


x 3 - r sinL . 


(47) 


Expanding Equation 46 and introducing the new constants and the transformation of Equation 47, 
makes Equation 46 become 

( F 2o)o ( 3x 3) 2 ~{( x i) 2+ ( x 2) 2+ ( x 3) 2 }] 2(F 21 ) 0 Xj x 3 2 (F 21 )j x 2 x 3 3( f 22)0 [( x i) 2 _ ( x 2 ) 2 ] 6(F 22 ) lXl x 2 2I 0 ~3I r 

2 ~ 7* + T 2 + r 2 + r 2 + r 2 ' 2 

(47a) 
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But, from their definitions (Equations 20b and 20c), 




[W 1 *(«.)’* («.)*] 

2I 0 = (A + B + C) ~ 2 


(47b) 


and 


[a(Xj ) 2 +b(x 2 ) 2 +c(x 3 ) 2 ] 

^ = 7 ? 

Equation 46 thus becomes (with Equations 47b and 47c in 47a) 


(47c) 


( X l) : 


( F 2 o) 0 


3 ( F 2 2 )o 


( X 2 ) : 


( F 2o)o 


3 ( F 2 2 ) 0 


K ) 2 ( F 


20) 0 +X l X 2 


[ 6 ( F 2 2 )l] +X 1 X 3 [ 2 ( F 2l)o] 


+ X 2 X 3 [ 2 ( F 2 2 )l] 


( X l) 


A + B+C 
2 



A + B + C 
2 



A + B+C 
2 


3C 

2 


(48) 


Since all the terms in x ; x ; on either side of the equality in Equation 48 are linearly independ- 
ent, the only way Equation 48 can be satisfied identically for all sets of x 2 , x 2 , x 3 with non-zero 
coefficients is for 


( F 2o)o 

2 

, , (B + C-2A) 

+ 3(F 22 ) 0 2 

(49) 

( F 2 0 ) 0 
2 

_ ,/ F x - (A + C - 2B) 

3 \ F 2 2 ) 0 2 

(50) 

( F 

V (A + B-2C) 

20/0 ~ 2 

(51) 

0, which implies <p 2 2 ~ 0° or 180° , 

(52) 


2 ( F 2 i)o s 0 

(53) 


2(F 21 ), - 0 • 

(54) 


Thus, Equations 52, 53 and 54 determine (f 22 ), 3 (F 21 ) 0 = (f 21 ), - 0. 

Substituting (F 20 ) 0 from Equation 51 into Equation 49 determines (F 22 ) 0 as 


(^22) 0 



(55) 
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It may be verified that Equation 50 is satisfied identically with (F 22 ) 0 and (f 20 ) 0 given by Equa- 
tions 55 and 51. 

Therefore, with respect to (r, L, <t>), or with reference to the principal axes of the mass dis- 
tribution, the potential (Equation 43) becomes 

+ order (-f) + •*• • (55a) 

In terms of the coefficients of Reference 5 this may be written as 


r 


GM 


( F 2o)o 


2 — (3 sin 2 L - l) - 3(F 22 ) 0 cos 2 L cos 20 


GM e 

r 


1 - ( a / r ) 2 (3 sin 2 L _ l) - 3J 22 (a/r) 2 cos 2 L cos 2 0 


+ orde 




■ (56) 


The constant a will be defined later as a mean equatorial radius of the model earth ellipsoid. 
In Equation 56 (from Equations 51 and 55 in Equation 55a) 


and 


( F 2o)o = 2C - (A. + B) 
M 3 (I) 2 ' 2M e (a) 2 


(^ 22)0 _ A ~ B 

M e (a) 2 " 4M e (I) 2 ‘ 


(57) 


(58) 


Equations 57 and 58 thus provide links between the gross mass properties of the earth (its principal 
moments of inertia) and its external gravity field. In fact, the oblateness of the earth (measured by 
2C- [A + B] ) is known more accurately from satellite determinations of j w than from direct surface 
measures of oblateness by astronomical observations of the change of latitude with distance north. 
Similar links can be established between the higher order gross mass properties of the earth and 
the higher order coefficients of the external gravity field. In part VII these coefficients will be 
linked to the surface geometry of the geoid, or mean sea level of the earth. 

The gravitational equipotential of Equation 56 with the higher order terms neglected and a 
centrifugal term added will be shown in part VII to give essentially an ellipsoidal surface if j 20 
and J 22 are sufficiently small. j 20 and j 22 may then be determined as functions of the rotation, 
oblateness and ellipticity of a model triaxial- ellipsoidal earth whose surface would be at mean sea 
level. 

Two questions remain to be resolved before giving geometric (geoid) significance to J 20 and 
J 22 . The first concerns the validity of ignoring the l/V 2 term in the gravitational potential for 
those field points (i.e., on the surface of an ellipsoidal earth away from the equatorial major axis) 
for which some of the mass distribution is farther from the origin than the field point. The infinite 
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series of harmonic functions in Equation 41 is the general solution of Laplace’s equation with the 
boundary condition V e (r -*°°) - o and represents the potential at all points exterior to the surface 
of the mass distribution. (See Reference 1, Chapter 13, for a rigorous proof of this statement.) But 
at points more strongly exterior to the surface the direct integration by series for the potential 
shows that the l/r 2 term is absent for the choice of origin at the c.m. The numerator of the l/r 2 
term in the harmonic series of Equation 41 is a function of l and only. On a radial line, there- 
fore, it is constant. Since it is zero in the outer region, it is zero everywhere to the surface be- 
cause the inner region to the surface can be completely covered by radial lines from the enfolding 
outer region. 

The second question concerns the orientation 
of the x p x ^ , x 3 principal axes of inertia in the 
earth. Let us assume the earth is spun about 
an arbitrary axis in inertial space. From Ref- 
erence 5, pp. 377-378, it is shown that for a 
rigid body such an initial rotation is invariable 
in the absence of external torques if, and only 
if, that initial spin axis is a principal axis of 
inertia. Astronomical observation shows that 
to a high degree of accuracy the earth’s north 
pole spin axis is invariable in direction in 
inertial space after accounting for all external 
torques. We can thus arbitrarily (between x p 
x 2 and x 3 ) identify x 3 as the polar axis of ro- 
tation of the earth. In the model to which the 
simplified series in Equation 56 belongs, L cor- 
responds to geocentric latitude with respect to 
the c.m., and </> is the geocentric longitude measured counterclockwise (looking down on the north 
pole) from the principal axis on the equator (Figures 2 and 3). 


(NORTH POLE) 



VI. MATCHING THE LEVEL SURFACE OF THE SECOND-ORDER GRAVITY AND EARTH ROTATION 

POTENTIAL TO THE SURFACE OF AN ELLIPSOID OF THREE AXES (THE TRIAXIAL GEOID] 

In this section we will show that the mean sea level surface of the earth, the geoid, is, to second 
order, an ellipsoid of three axes. Arbitrarily the x t axis is taken along a, the semimajor equatorial 
axis of the triaxial ellipsoid. 

In the ellipsoidal earth model of Figure 3, e is the ellipticity of the equator and f is the polar 
flattening; and in terms of a , b , and c , the half axes of the ellipsoid, 
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and 


where 


f = l 


c a - c 


a + b 


(60) 


(61) 


^4iiaiiuuS uu olid 61 oT6 arbitrary but convenient definitions Of 1 and a.) The 5 Ui7 1 a C e of the el- 
lipsoid in Figure 3 is given by 

(t) + (tt) + (t) = 1 (62) 

In terms of the coordinate transformation of Equation 47 the surface is also given by 

x = R cos L cos 0 


R cos L sin 0 


R sin L 


(63) 


where R is the radius to the surface of the ellipsoid. 
Combining Equations 62 and 63 produces 


cos 2 L cos 2 0 cos 2 L sin 2 0 sin 2 L , „ 

* Si * — = Vr» . 


(64) 


Multiplying Equation 64 through by b 2 , solving for b/a from Equation 59 and rearranging makes 
Equation 64 become 

( 1 - e) 2 (cos 2 L cos 2 4 >) + cos 2 L sin 2 0 ~ b 2 . (65) 

Next solve for c 2 from Equation 60 and introduce this into Equation 65 which then becomes 


b 2 

( 1 - e) 2 cos 2 L cos 2 0 + cos 2 L sin 2 0 = rrr 


R 2 


sin 2 L 
(1-f) 2 


( 66 ) 


A combination of Equations 59 and 61 leads to 


bf _ [2(1 -e)1 

a 2 " L 2-e J 


( 67 ) 
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Upon the substitution of Equation 67 into 66 this equation now becomes 


(i)' 


( 1 - e ) 2 cos 2 Lcos 20 + cos 2 L sin 2 0 + S * n \ 

( 2 -e)(l-f ) 2 


" 2 ( 1 -e) : 
2~e 


(68) 


When Equation 68 is somewhat simplified, it takes the form 


( ^/R) 2 - 4 ; ( 2 - e ) 2 cos 2 L cos 2 0 + 


4 (2 ~ e ) 2 cos 2 L sin 2 0 s in 2 L 


( 1 ~ e ) 2 


(1 - f ) 5 


(69) 


Two factors from terms in Equation 69 are now expanded to give 


2 -e 


1 - e 


® ( 4 - 4e + e 2 ) (l + 2e + 3e 2 + • * •) - 4 + 4e + 5e 2 , for e << 1 


(70) 


and 


( 1 - f )“ 2 - 1 + 2f + 3f 2 , for f«l 


(71) 


With the approximations of Equations 70 and 71, Equation 69 becomes 


(a/R) 2 - 4(2~e) 2 cos 2 Lcos 2 0+^;(4+4e+5e 2 ) cos 2 L sin 2 cp + (l + 2f + 3f 2 ) sin 2 L • (72) 

Equation 72 describes the surface of the triaxial ellipsoid to e 2 and f 2 in the equatorial ellipticity 
and oblateness (polar flattening). On the surface of the rotating earth ellipsoid the centrifugal force 
per unit mass is 


(fJr - ( V r') R f ' = (" 2r ') R f' • (73) 

where v is the surface velocity at r' from the axis of rotation of the earth, and o) is the rotation 
rate of the earth. The force field of Equation 73 is derivable as the gradient (V = r' d/dr*) of a 
potential function of r': 


( V c) R 


(^ r ') 2 R 
2 


(74) 
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The centrifugal potential at the surface of the ellipsoid is thus given as 
2 2 

[ V c]r = ^2 [( x i) 2 + ( x 2) 2 ] r ~ ^2 [ r2 cos 2 L cos 2 (p i R 2 cos 2 L sin 2 0] , 


~2 R 2 cos 2 L 


( 75 ) 


By combining Equations 56 and 75, the full potential at the surface of the ellipsoid can be ex- 
pressed as 


V( total ) 
where 


= v +v„ 


= ™ jf“ + \ ^ ) cos 2 L- — Y (§■) (3 Sin 2 L - l) - 3J 22 cos 2 Lcos2^j> (76) 


QA 


(77) 


Next, the assumption is made that the surface of the earth ellipsoid is a level potential surface 
of constant amount c 0 , where (from Equation 76) 


t[t + !(D cos^L- S -f (§) (3sin 2 L-l) -3J 22 (r) 


cos 2 L cos 20 


(78) 


Solving for R from Equation 72 will give 


R = a 


jl + e + cos 2 L cos 2 0 + ^ ^ cos 2 L sin 2 <p + (2f + 3f 2 ) 


sin 2 L 


- 1/2 


which can be simplified by neglecting all terms in e and f above first order. This simplification 
leads to 


R = a [l - e cos 2 L cos 2^ + 2f sin 2 l] 1/2 . (79) 

The following expressions are then obtained by an expansion of Equation 79 (to first order in f and e) : 


R" 1 

R 2 


a ^1 + cos 2 L cos 20 - f sin 2 LJ , 

(79a) 

(a) -1 ^1 - cos 2 L cos 20 + f sin 2 , 

(80) 

(a) 2 (l + e cos 2 L cos 20~2fsin 2 L) , 

(81) 
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and 


R 3 


(a ) 3 


M 


cos 2 L cos 20 - 3f sin 2 L 


(82) 


When Equations 80 through 82 are substituted into Equation 78 and the second-order product terms 
in J 20 , J 22 , e, f and a are ignored, Equation 78 takes the form 


C„ = — 


GM 

a 


cos 2 L cos 20 + f sin 2 L + tt 


(l - sin 2 l) - (3 sin 2 L - l) “ 3J 22 cos 2 L cos 20 


? . 


(83) 


In order for Equation 83 to be an identity for all L and 0 on the surface of the earth ellipsoid, the 
following conditions must be satisfied: 


C 0 



a 

2 + 



(84) 


2 3J 22 - 0 , 


(85) 


_ a _ 3 

' 2 2 J20 


0 


( 86 ) 


The existence of the geometric-potential compatibility Equations 84-86 establishes 
of the surface equipotential of second order gravity and earth rotation as a triaxial 
Equation 85 it is seen that 

J 22 = -e/6, 


approximation 
ellipsoid. From 

(87) 


and from Equation 86 the gravity potential constant J 20 is identified as 


J 


20 


2f _ *_ 

3 “ 3 ' 


( 88 ) 


Expressions 87 and 88 give approximate physical interpretation to the harmonic coefficients in the 
potential of Equation 56 as functions of the small parameters of the model earth triaxial ellipsoid 
and its rotation. 

This now allows us to describe the geometry of the earth model and identify axes x r and x 2 from 
satellite and surface derived values of J 20 and J 22 . 

According to Kozai (Reference 7) and Wagner (References 8 and 11) the best values of J 20 , J 22 , 
GM and a are 


J 20 

= 1.08248 x 10" 3 (Kozai) , 

(89) 

J 22 

= - 1.9 x 10~ 6 (Wagner) . 

(90) 
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From O’Keefe, Eckels and Squires cited in Reference 5, the following constants are obtained: 


GM - 3.98603 x 10 s km 3 /sec 2 (Kozai) 

a = 6.378165 x 10 3 km 

co = -7292115 x 10" 4 rad/sec (from Reference 9) 

Therefore, from Equations 91 and 92 in Equation 77 


(91) 

(92) 


3.46165 x 10‘ 


(92a) 


Solving for the oblateness (polar flattening) of the ellipsoid model from Equation 86 with the con- 
stants of Equations 89 and 92a gives 


f * 1.5 x 1.08248 x 10~ 3 + .5 x 3.46165 * 10“ 3 « 1/298.24. (93) 

Similarly, a solution for the equatorial ellipticity from Equation 85 with j 22 given in Equation 90 
will give 


e = 11.4 X 10“ 6 . (94) 

The difference in major and minor radii of the elliptical equator of the model triaxial earth can 
now be found. From Equations 59 and 60 


a( 1 - e) - b - 0, 

a + b = 2a • 

Solving Equations 95a and 95b for a and b will show that 


(95a) 

(95b) 


a - 2a/ ( 2 - e ) , 


2a(l - e) 
2 - 


Therefore, from Equation 96 


(96) 


2a e ae 

a " b ” T^T ~ 1 - e/2 


ae 


(97) 


With the results of Equations 91 and 94 in Equation 97 the difference in major and minor radii of 
the elliptical equator is calculated as 235 feet. 
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The order of the principal axes of inertia x 1? x 2 and x 3 is established by the signs of j 20 and 
j 22 in addition to their relative magnitudes. Since J 22 is negative, from Equation 58 it can be seen 
that B > A. Also from Equation 57, since J 20 is positive, 2C > (A + B) . Therefore since B > A, 2C > 2A, 
or C> A. The axis of minimum moment of inertia, x 1? is thus established. On the earth ellipsoid 
(Figure 3) this is the major axis of the equatorial ellipse. (Note also for j 22 <0, e>o from Equa- 
tion 85 and therefore a > b from Equation 59.) To establish the probable order of principal axes x 2 
and x 3 , consider the earth ellipsoid as homogeneous. Then (from Reference 6, p. 292) 


(c 2 /a 2 ) + 1 

,/c = (SM7- 


so that if we can show b/c > 1 , then C > B. But from Equation 60 c = a( l - f ) . This result with that 
from Equation 96 will give 


b 

c 


(i- f)d-f> 



(from Equations 93 and 94). This final result establishes the absolute order of the moments of 
inertia (for a homogeneous ellipsoid) as 


A < B < C , (98) 

and for the actual earth, as least A < B , c . 

VII. THE LEVEL SURFACES OF THE GRAVITY POTENTIAL 
OF THE EARTH TO FOURTH ORDER 

The form of the potential function in Equation 56 (and the magnitude of the second-order coef- 
ficients) suggests the interpretation that the earth attracts principally as a sphere (homogeneous 
in spherical shells). The higher order terms may be interpreted as arising from small mass de- 
viations in the principal spherical model. At the surface these deviations are reflected as slight 
corrugations in latitude, longitude and radius above and below the surface of the average earth 
sphere. In particular, the previous development has shown how the principal term together with 
the j 20 and J 22 harmonics and the centrifugal potential describe an earth with an equipotential 
surface as a triaxial ellipsoid only very slightly deviating from a spherical shape. In the following 
discussion the surface deviations from an average earth sphere of radius a (chosen for convenience) 
due to the higher harmonics will be considered. The deviations will be such that they give an 
equipotential surface for the particular harmonic. The actual equipotential surface of the earth is 
then composed of the higher harmonic deviations superimposed on the second order triaxial el- 
liposid described in the previous section. Higher harmonics to the fourth order will be considered. 

From Equation 42, since 4>~4> n ° \-k , f. A - F. . n 0 and F 91 - 0 as discussed in Sections V 

nm nm "10 11 ^ A 

and IX, the geopotential through the j 44 term with respect to geocentric-geographic coordinates 
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r^-L, X. measured east from Greenwich (Figures 2 and 4) is: 


1 - (i/r) 2 J 20 P 2 ° (sinL) - (a/r) 2 J 22 P 2 2 (sinL) cos 2 ~^ 22 ) 

- (a/r) 3 J 30 P 3 ° (sinL) - (i/r) 3 J 3l P j (sinL) cos (^~^ 31 ) 

- (a/ r ) 3 J 32 P 3 2 (sin L) cos 2 (X - \ 32 ) - (a/r ) 3 J 33 P 3 3 (sin 1 ) cos 3(^ - A- 33 ) 

- (a/r) 4 J 40 P 4 ° (sinL) - (a/r) 4 J 4 , P ^ (sinL) cos (X-A- 41 ) 

-(a/r) 4 J 42 P 4 2 (sinL) cos 2(K - _ (i/r) 4 p^3 ( s inL) cos 3(\ - \ 43 ) 

" (a/ r ) 4 J 44 P 4 4 (sin L) cos 4(\ - ^ 44 )J ■ (99) 

where 


V e (geocentric -geographical ) 


CM 

r 


P 


m 

n 


"d" P n (sinL)' 

~ cos'” L 

_ d(sin L) m 


The K above represents the geographical longitude of the principal plane of longitudinal symmetry 


for the nm harmonic. 

Equation 99 is Equation 42 with the coef- 
ficients rewritten to conform with the J coef- 
ficients of Equation 56. The P n m are the H nm of 
Equation 36 without the coefficient D nm . With the 
selection of the reference coordinate system as 
explained in Section V, the l/r 2 term in the 
potential and the J 21 term have been ignored. 
Depending on the actual choice of geocentric- 
geographical coordinates of radius and latitude 
used in a given application, reported coefficients 
for Equation 99 may vary slightly. For example 
if a calculation is to be run with Equation 99 in 
geographic coordinates referenced to the north 
pole and an earth center on the equatorial plane, 
the reported coefficients would appear to be 
logically in order if the model used to derive 
those coefficients used the same references. 
However, since the north pole experiences small 
anomalous motions which are generally not 


PROJECTION ON 
EQUATORIAL PLANE OF 
P (TEST MASS) 



NOTES : 

$ “ fnm = X “ X nfn 


Figure 4 — Longitude relationships in the gravity potential 
of the earth with respect to the elliptical equator of the 
earth. 
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corrected for in the calculations, very small discrepancies may still logically exist between other- 
wise correct sets of coefficients. More basic still is the consideration that the uncorrected polar 
spin axis is not quite an axis of principal moment of inertia for the earth. For calculations using 
an equatorial center and latitude with respect to the spin pole, a "best" set of coefficients should 
include the J 21 term (see Section IX). 

In what follows the potential of Equation 99 is referred to an origin at the c.m. of the earth 
and the north pole spin axis, which is assumed to be a principal axis of inertia for the earth. 

Table 1 gives longitude coefficients (j nm , k nn9 m / 0) for the expansion Equation 99 as reported 
by geodesists from 1915 to 1964. 

The following material in this section describes how this so-called spectral representation of 
the field can be translated into a spatial representation of deviations on the mean earth geoid. 
Summing up each harmonic spatial deviation of the geoid finally gives the actual sea level surface 
described by that particular set of harmonics. Figures 16 through 23 give the actual sea level 
surfaces for some of the more recent longitude geoids in Table 1. The zonal harmonics (j 20 , J 30 > 
J 40 ) have not been reported in this table for each investigation since they do not influence deviations 
of the geoid along any given latitude circle. Figures 16 through 23 show that the full surface 
equipotential contours are sensitive to these harmonics however. The reference geoid in each case 
has a flattening of 1/298.2 which includes the deviatory effect of the centrifugal potential and the 
oblate gravity term J 20 (Section VI). Though not listed, all of these recent geoids (Figures 16 through 
23) report J 20 , J 30 and j 40 values little different from the set given by Y. Kozai in Table 1 (also 
Reference 7): 


J 20 = 1.08248 x 10“ 3 

J 30 = -2.6 x 10“ 6 

J 40 = -1.84 x 10“ 6 . 

A truer comparison of the geoids in Table 1 can be gained through the integrated spatial repre- 
sentation such as Figures 16 through 23. This is because it is the full field at various locations that 
is actually measured in any investigation. This representation illustrates the fact that relatively 
minor spatial differences can give rise to large spectral differences term-by-term. It also points 
up the necessity of considering the reported spectral coefficients as a set in any application, at least 
for the relatively ill-determined longitude ones. Considered as a set, or spatially, it will be seen 
that most of the recent geoids agree rather well in overall features. In theory the spectral effects 
are well separable (as in practice they have proved to be for the zonal terms). In the longitude 
reductions to date, the data has been too lacking in quantity and quality to allow meaningful separa- 
tion except perhaps by order of magnitude. The one exception to this state of affairs is in regard 
to J 22 (see Reference 11, for example). Nevertheless, on the basis of the spatial representations 
of Figures 16 through 23 derived from different kinds of gravity data, it is possible to delineate a 
number of large area features over the earth (sensitive to the longitude harmonics) which appear in 
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Longitude Coefficienta in the Gravity Potential of the Earth i V E = -p 2 _] r~ ( R o/ r ) " p n"" ( sin <£) J nm cos m (\ - \ nm ) 11 , as Reported 1915-1964*. 
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cemmon. One or more of these common longitude-sensitive undulations of the geoid can be said to 
exist with the same certainty as the mass deviations responsible for the so called "pear-shaped 
earth." 

Table 2 lists ten such common features seen in the spatial representations of Figures 16 
through 23. It is noted that only two appear to be absolute in the sense of standing out strongly on 
all the geoids (see Table 2A). The reader is welcome to try his own eye for comparative land- 
scaping in these figures. Perhaps the feature called "the middle American low" ought to be com- 
bined with "the central Pacific and central Atlantic low" features to distinguish a more common 
strong tendency of these geoids. The same might be said of "the Afro- European" and "the Afro- 
Antartic" highs. 


Table 2 

Longitude Geoid Features (see Figures 18-23). 


Geoid (Figures 16-23) 
(References are in 
Table 1 except as 
noted.) 

Eastern Hemisphere 

Western Hemisphere 

Afro- 

Antartic 

High 

Afro- 

European 

High 

Indian 

Ocean 

Low 

Micronesian 

High 

Austro- 

Antartic 

Low 

Alaskan 

High 

Central 

Pacific 

Low 

South 

American 

High 

Middle 

American 

Low 

Centra! 

Atlantic 

Low 

Izsak (1964)* 

/ 

p 

/ 

/ 

/ 

p 

/ 

✓ 

p 

/ 

Guier (1963) 

/ 

/ 

/ 

/ 

✓ 

p 


? 



Kaula (Sept. 1963) 

/ 

/ 

/ 

/ 

✓ 


p 

p 

p 


Izsak (1963) 

/ 

✓ 

✓ 

/ 

✓ 


p 


p 

p 

Anderle and Oester- 
winter (1963) 

/ 

/ 

i/ 

/ 

p 

/ 

p 

✓ 

/ 


Kozai (1962) 

/ 

✓ 

✓ 

/ 

✓ 


✓ 

✓ 

p 

✓ 

Uotila (1962) 

p 

✓ 

✓ 

/ 

/ 

p 

✓ 


✓ 

p 

Kaula (June 1961) 

/ 

/ 

/ 

✓ 

✓ 

p 

p 

✓ 

p 



•See Reference 10. 


Table 2A 

Longitude Geoid Features in Order of Common Strength. 

1. Indian Ocean Low 

2. Micronesian High 

3. Afro-Antartic and Afro- European High 

4. Austro- Antartic Low 

5. South American High 

6. Central Pacific - Central American - Central Atlantic Low 
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The deviations giving rise to the terms of Equation 99 are described below. 


a. Deviations of the Spherical Equipotential Surface Associated with the J 2 o Term of the 
Gravity Field of the Earth 

To consider the effect on the average earth sphere (of radius about equal to a) of just the J 20 
term in Equation 99, we write the potential of that sphere modified by the j 20 deviation as 

V - ~ jl - (a/r) 2 — ' 2 ~ (3sin 2 L~l)| • (99a) 

According to Y. Kozai in Reference 7, J 20 = 1082.48 x 1 0 6 . The deviations from the average 
sphere are referenced to latitude L = ±35.2° since the deviatory part of Equation 99a goes to zero 
there. Thus when L = 35.2° in Equation 99a, r = a. The equipotential surface for Equation 99a 
then has a potential of constant amount GM/a. The condition equation for the equipotential surface 
due to Equation 99a is thus: 


CM 

a 



^ J 2 o( 3sinJL “ x ) ' 


(99b) 


Let 


r - a + Ar 20 


(99c) 


where 


Ar 


20 


« a . 


Next introduce a + Ar 20 into Equation 99b and solve for Ar 20 to get 


Ar 


20 


J2O 

2 


! L “ l) 


The zeros of Equation 99d are at L - ±35.2° (see Figure 5). Extreme radial deviations due to 
are attained when 


(99d) 

J20 


^20 

dL 


0 - - 3a S * 2 ^ ’ or when L ~ 0° and ±90° 


The extreme radial deviations at these latitudes are 

( Ar 2o) min = " aJ 2 o ( at the poles) 

= - 6378 x 1082 x 10" 6 = -6.90 km. = -4.28 st. miles ; 
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void 


( Ar 2o) ma , 


3 J20 

2 (at the equator) 


= 2.14 st. miles (Figure 5). 


The deviations in Figure 5 are strictly 
gravitational. When the surface rotational po- 
tential of the earth is considered as well, an 
equipotential surface with the approximate shape 
of a biaxial ellipsoid is the result (ignoring 
much smaller higher order effects) as Section VII 
showed. 



Figure 5— J 20 deviations from an average earth sphere 
from a potential from Reference 7. 


From Equation 60 the polar radius c of the oblate earth spheroid is given in terms of 
the flattening f and the mean equatorial radius I as 

c = a(l - f ) . 

The difference, then, in polar and equatorial radii for the oblate earth spheroid (often used as the 
ellipsoid of reference for the earth or the "geoid") is 

a - c = a( 1 + f - 1 ) = af = 6378/298.2 (from Equation 93) = 21.4km . 

Thus, the maximum deviations of the oblate earth spheroid from the average earth sphere are of 
the order of 21.4/2 = 10.7 km = 6.65 st. miles when the rotational potential is considered as well 
as the deviatory part of the gravitational potential. This latter result is consistent with the maximum 
and minimum equipotential deviations calculated above for just the gravitational part (due to J 20 ). 

6. Deviations Associated with the J 22 Term of Equation 99 

(According to C. A. Wagner in Reference 11, J 22 = -1.9 x 10' 6 , ^ 22 = -20°.) The J 22 deviatory 
geopotential term in Equation 99 is 

(V e ) 2 2 = -3(j 22 )(I/r) 2 cos 2 Lcos2(\-^ 22 ) • (99e) 

Proceeding as in the case for the deviations giving rise to the J 20 term, we arrive at radial devia- 
tions from the surface of the average earth sphere of radius a given by 


Ar 22 = -3a J 22 cos 2 Lcos 2(X - \ 22 ) . (99f) 

The zeros of Equation 99f are at L = ±90°, and at every latitude where a. = \ 22 ±45°, \ 22 ± 135° 
(Figure 6). 
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Extreme deviations due to J 22 are reached 
when 


d&r 22 

- 0 = 3a J 22 sin 2Lcos 2(X - X 22 ) 


and, 


aAr 2: 


r 0 - 6a J 2 2 cos 2 L sin 2(X - X 22 ] 


Figure 6-J 22 deviations of the mean earth sphere for simu i t aneously. These minimax deviations 
J 22 <0, from Reference 11. J 

occur (other than at the poles which are inflec- 
tion points where Ar 22 = 0) at the equator at longitudes where x - x 22 , X 22 ± 90°, x 22 + 180°. The 
absolute magnitude of these extreme deviations is 


|Ar 22 (max. )| = -3aJ 22 = 3 x 6378 x 1.9 x 10~ 6 = 36.4 meters 

= 119 feet . 


This corresponds to a difference in major and minor equatorial radii of 238 feet. Compare this 
with Equation 97 (also see Figure 6). 


c. The Deviations of the Surface of the Earth Associated with the 
J 30 Term of Equation 99 

The effect on the average earth sphere, whose radius is about equal to a, of just the j 30 term 
can be calculated when the potential of that sphere is modified by the j 30 deviation to give 


v 


GM 

r 


1 " (a/r) 3 



(5 sin 3 L ” 3 sin L ) 


( 100 ) 


(According to Y. Kozai in Reference 7, J 30 = -2.6 x 1(T 6 .) The deviations from the average sphere 
will be referenced to the equator since the deviatory part of Equation 100 goes to zero there. Thus, 
when L = 0, r " a in Equation 100. The equipotential surface for Equation 100 thus has constant 
amount o^I. The condition equation for this equipotential surface is then 


GM 

a 


GM 

r 1 


(a/r ) ; 


¥ (5 


in 3 L ~ 3 sin L ) 


( 101 ) 
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When Equation 101 is rearranged to solve for r, it becomes 


(r/a) 4 - (r/a) 3 


( 5 sin 3 L - 3 sin L ) 


( 102 ) 


Let r = a + Ar^, where Ar 30 « a . Then, to the first order in Ar/a, the left side of Equation 102 be- 
comes Ar/a and Equation 102 is written approximately as 


Ar 


30 



($ sin 3 L - 


3 sin l) . 


(103) 


The zeros of Equation 103 (besides L = 0) are at L = ±50.7°. Furthermore, Ar 30 < 0, 01Li50.7°and 
-50.7° < L < -90°; Ar 30 > 0 everywhere else (see Figure 7). 


Extreme radial deviations are attained when 


d ( Ar 3o) 

dL 


a h 


( 15 sin 2 L cos L _ 3 cos L ) - 0 


(104) 


The zeros of Equation 104 (besides L = 0) are at L = ±26.6° and at 
L = ±90°. 



-50.7° 


(MIN. 54') 


When the value of a in Equation 91 and J 30 as quoted above are used, 
Ar extreme, calculated from Equation 103, is equal to 24 feet at l = ±26.6° 
and 54 feet at the poles (Figure 7). 


Figure 7— J 30 deviations 
(Reference 7). 


d. Deviations Associated with the J 31 Term of Equation 99 
The J 31 deviatory geopotential term is 


“J 31 (a/r) 3 -j (l5 sin 2 L ~ 3) cos L cos (^ ” ^ 31 ) • 


By proceeding as in the case for the deviations giving rise to the J 22 term, the radial deviations 
from the surface of the average sphere of radius a can be shown to be 


^r 31 - -a ~ ?r (l5 sin 2 L - 3) cos Lcos ~ ^ 31 ) 


(105) 


The zeros of Equation 105 are at 

L = ±90° (all longitudes) 

L - ±26.6° (all longitudes) 

K - 90° + ^ 31 and 270° + ^ 31 (all latitudes) (Figure 8). 
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Extreme deviations due to J 31 are reached when 



Figure 8— J 31 deviations for J 31 < 0 (Appendix A). d K 


(l5sin 2 L-3)sinL 
+ 30 sin L cos L J " 0 ,(106a) 


- sin(x - \ 31 )(l5 sin 2 L ~ 3)cos L - 0. 

(106b) 


The simultaneous solution of Equations 106a and 106b yield two sets of extreme deviation points 
across the average sphere. One is at the crossing of the "nodal lines" of Figure 8 where Ar is 
zero. These are evidently flat points of inflection for the equipotential surface due to J 31 . The 
other set is at the intersection of \ - \ 31 and 180° + \ 31 with L = ±70° and 0°. Four of the six 
extreme points of this set yield the same absolute deviation of -1.75 a J 31 . The other two yield 
identical absolute deviations of -3 a j 31 /2 (Figure 8). 


e . Deviations Associated with the J 32 term of Equation 99 
The J 32 deviatory geopotential term is 


-J 32 (a/r) 3 15 cos 2 L sin L cos 2(^ _ 2 ) .• 

Deviations are calculated as in d. above, and 

Ar 32 - ~a J 32 15 cos 2 L s in L cos _ Si) (107) 

with the results appearing on Figure 9 below. 


/. Deviations Associated with the J 33 Term of Equation 99 
The J 33 deviatory geopotential term is 


EAST 


-J 33 (a/r ) 3 15 cos 3 L cos 3(A. - \ 33 ) 

Deviations are calculated as in b . above and 

Ar 3 3 "a J 33 15 cos 3 L cos 3(^ ' ^ 33 ) (108) 
with the results appearing in Figure 10 below. 


Ar 32 (MIN. ) =2.885 J 32 



Ar 32 (MAX. ) = “2 .885 J 32 

Figure 9— J 32 deviations for J 32 < 0 (Appendix A). 
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g. Deviations Associated with the 
J u0 Term of Equation 99 

The J 40 geopotential term is 

— J 4 o (5/r) 4 [^T sin 4 L- 30 sin 2 L + 3 ] . 

(According to Reference 7, j 40 = -1.84 x 10" 6 .) 
Deviations are calculated as in a. above and 


Ar(MIN. ) = I 5 SJ 33 



Figure 10— J 33 deviations for J 33 < 0. 


Ar 40 - "a J 40 (y sin 4 L - 30 sin 2 L + 3^(109) 
with the results appearing on Figure 11. 

h. Deviations Associated with the 
J l)] Term of Equation 99 

The J 41 deviatory geopotential term is 

■J 41 (a/r) 4 ('§')( 140 sin 3 L 

- 60 sin l) cos L cos ( X. - X 4 ^ . 
Deviations are calculated as in d. above, and 

Ar 4 i = ~a J 41 (-g-)(l40sin 3 L 

- 60 sin l) cos L cos ( X - X 41 )• (110) 

with the results appearing in Figure 12. 


(MAX. 38') 


\ 



/ 


(MAX. 38') 

Figure 11 — J 40 deviations 
(Reference 7). 


Ar 41 (ABSOLUTE MAX. )= -2.63a J 41 



Ar 41 (RELATIVE MIN. ) = 1 .81 a J 41 

Figure 1 2 — J 4 1 deviations for J 41 < 0 (Appendix A). 


i. Deviations Associated xvith the J u2 Term of Equation 99 
The J 42 deviatory geopotential term is 


-J 42 (a/ r ) 4 (^r)( 420 sin2 L “ 60) cos 2 L cos 2(X - ^ 42 ) 
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Ar 42 (ABSOLUTE MAX. ) = -9.65 J 42 



Ar 42 (RELATIVE MIN . ) = 7 .5a J 42 
Figure 1 3 — J 4 2 deviations for J 42 < 0 (Appendix A). 


Ar^ (MAX. ) =-34.15J 43 



Ar^ (MIN.) = 34. 1 a 

Figure 14 — J 43 deviations for J 43 < 0 (Appendix A). 



Ar^ (MIN. ) = 105a 


Figure 1 5 — J 4 4 deviations for J 44 < 0 (Appendix A). 


Deviations are calculated as in d. above, and 

Ar 42 = -Ij 42 (4-)(420 sin 2 L 

- 60 ) cos 2 L cos 2 (A - \ 2 ) (in) 

with the results appearing here on Figure 13. 

j. Deviations Associated with the 
J. , Term of Equation 99 

The J 43 deviatory geopotential term is 

-J 43 (a/r ) 4 (~g-) 840sinLcos 3 Lcos 3(K-^ 43 ). 

Deviations are calculated as in d. above and 

Ar 43 = -a J 43 (4") 840 sinLcos 3 L cos 3(A - \ 3 ) 

( 112 ) 

with the results appearing here on Figure 14. 

k. Deviations Associated with the 
J !f!j Term of Equation 99 

The J 44 deviatory geopotential term is 

-J 44 (a/r) 4 (-^) 840 cos 4 L cos 4 (A -\ 4 ) 
Deviations are calculated as in b. above and 


Ar 44 = -a J 44 (-§-) 840cos 4 Lcos4(A-\ 4 ) (113) 

with the results appearing on Figure 15 above. 

The complete geocentric-geographic gravity potential of Equation 99 with the zonal coefficients 
as reported in Reference 7 and J 22 , \ 22 as given in Reference 11 is given below. The longitude (\) 
reference is Greenwich with positive angles to the east. The radial (0 reference is the center of 
mass of the model earth. The latitude (l) reference is the equatorial plane of the model earth 
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parsing through the c.m. and the north pole spin axis. Thus the spherical harmonic representation 
of the gravity potential of the earth through fourth order takes the form of 

v e = “7“ 1^1 - ( 1082.48 x 10~ 6 )(a/ r ) 2 ( 2 )(3 sin 2 L~ l) - (~1.9x 10~ 6 )(a/r) 2 3cos 2 L cos 2(k + 20°) 

- (-2. 56 x 10~ 6 ) (a/r ) 3 (^(s sin 3 L - 3 sinL)- (j 31 ) (a/r) 3 (-j) cos L (l5 sin 2 L~ 3 ) cos (\ ~ 

~ ( J 32 ) (a/r ) 3 15 cos 2 L sin L cos 2(\ - ^ 32 ) “( 133 ) (a/r ) 3 15 cos 3 L cos 3(\ - ^ 33 ) 

- (-1 . 84 x 10” 6 ) (a/r ) 4 (“g^lsin 4 L-30sin 2 L+3)~(j 41 ) (a/r ) 4 (~g“)(l40 sin 3 L~ 60sinL) cos ” ^ 41 ) 

_ (J 42 ) (a/ r ) 4 (s) (420 sin 2 L ~ 60 ) cos 2 L cos 2 ~ 2 ) - (j 43 ) (a/r ) 4 840 sin L cos 3 L cos 3 - ^ 43 

- ( J44 ) & X > 4 (f) 840 C ° S4 L COS 4 ( K ~ X 44 ) ] • ( 114 ) 


with a = 6378.165 km. 


GM_ = 3.98603 x 10 5 kmVsec 2 



Figure 16 — Sa tel life -Camera geoid from Izsak, 1964 (Reference 10) (Level curves of geoid heights 
at 10 m intervals with respect to a reference geoid of flattening 1/298.2), 
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Figure 23— Astro-geodetic ond gravimetric geoid from Kaula, June 1961 (Table 1) (geoid heights in meters). 

VIII. THE GRAVITY FORCE FIELD OF THE EARTH THROUGH FOURTH ORDER 

To display the gravity force field of the earth in geocentric coordinates referred to the earth's 
spin axis, we follow the general prescription that with respect to any orthogonal coordinate system 

F(gravity) = VV e (gravity) (115) 

(see Section n and References 1 and 2). 

The results of taking the gradient of v e from Equation 99 gives the force field of the earth 
through fourth order as 


F = fF + AF, + LF, = w 

r A L e 


dr r cos L ~dk + r dL 


_ f __ r 3 

F r = + (a/r ) 2 [*2 J 20 (3 sin 2 L- l) + 9J 22 cos 2 L cos 2(\ - \ 22 ) 

+ 2(a/r) J 30 (5 sin 2 L - 3 ) (sinL) + 6(R 0 /r) J 31 (5 s in 2 L - 1 ) cos L cos [k ”^ 31 ) 

+ 60(a/ r ) J 32 cos 2 L sin L cos 2(X. - \ 32 ) + 60 (R 0 /r ) J 33 cos 3 L cos - ^ 33 ) 

5 25 

+ 8 (a/r) 2 J 40 (35 sin 4 L ~ 30 sin 2 L + 3 ) + (a/r) 2 J 41 (7 sin 2 L “ 3 ) cos L sin L cos ( \ - \ 41 ) 

75 

+ ~~ 2 ~ (a/r ) 2 J 42 (7 sin 2 L - l) cos 2 L cos - \ 42 ) + 525(a/r ) 2 J 43 cos 3 L sin L cos 3(\ - \ 43 ^ 

+ 525(a/r) 2 J.. cos 4 Lcos 4(\ - \^)l\ ; (11 


F x = (a/r) 2 ^ 6J 2 2 cos L sin 2(\ - \ 22 ) +5 (a/r ) J 31 [5 sin 2 L - l] sin ( \ - - 31 ) 

+ 30(a/r ) J 3 2 cos L sin L sin 2 _ ^ 32 ) + 45(a/r ) J 3 3 cos 2 L sin 3(X - ^ 3 3 ) 

+ -^(a/r) 2 J 4I [7 sin 2 0“3]sinLsin[A.~X 41 ) + 15(a/r) 2 J 42 (7 sin 2 L - l) cos<£sin2(\~/V 42 ) 


+ 315(a/r) 2 J 43 cos 2 LsinLsin 3^X. - ^43) 

+ 420(a/r ) 2 J 44 cos 3 L s in 4(^ _ >V 44 ) ^ * (118) 

M e _ f , 

F l - — - (a/r ) 2 v '3J 20 sin L cos L + 6J 2 2 cos L sin L cos 2(^ - -V 22 ) 

3 3 

- ~2 (a/r ) J 30 (5sin 2 L'l) cos L + ^j (a/r)J 31 (15 sin 2 L - ll) sinL cos (A. “^31) 

+ 15(a/r)J 32 (3sin 2 L“l) cos L cos 2(^ _ ^ 3 2 ) 

+ 45(a/r ) J 33 cos 2 L sin L cos 3(^ - ^ 33 ) ~ ( a /r) 2 J 40 (7sin 2 L-3)sinLcosL 
+ \ (a/r) 2 J 4 , (28 sin 4 L - 27 sin 2 L + 3) cos (\ -X 41 ) 

+ 30(a/r) 2 J 42 (7 sin 2 L~4) cos L s in L cos 2 [K - \ 4 2 ) 


+ X 05 ( a/r ) 2 J 43 (4 sin 2 L - 1 ) cos 2 L cos 3(X - ^ 43 ) 
+ 420(a/r) 2 J 44 cos 3 L sin Lcos 4(\ - ^44)} ■ 


(119) 


The earth gravity force field given in Equations 117 - 119 and the potential field given in Equa- 
tion 99 are with respect to geocentric spherical coordinates fixed in the earth, referenced to its 
spin axis (the north pole) and its center of mass (Figure 24). 

It is instructive to note the relative influence of the orders of earth tesseral gravity at different 
altitudes above the earth 1 s surface. In general the higher one goes, the more the earth behaves as 
a point mass with a l/r 2 field. Physically, this asymptotic behaviour must occur because at in- 
creasing distances from the earth, the distances to each point in the earth from a fixed outside 
point become more nearly the same. The field of a single mass point in the earth is, of course, 
just a simple l/r 2 field. But in addition to this total behaviour since each order of gravity, mathe- 
matically, has a scaling factor of (a/r) with respect to the next lowest order; this means that as the 
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NORTH 


distance (r) from the earth increases the lower 
orders of gravity play an increasingly stronger 
role. 

This might be anticipated from the fact that 
the lower orders of gravity measure mass 
anomalies within the earth which have wider 

QcpdJ. diiuii, t ct otu^ uiouuiv/cS the SuuSS 

anomalies distributed with close separation, 
are at more nearly the same distance from the 
outside field point than the more widely sepa- 
rated anomalies of lower order. 



Figure 24 — Earth-fixed geocentric coordinate system 
reference for the earth’s gravity field. 


Thus, these higher order gravity anomalies should have a less and less influential role in the 
total nonhomogeneous earth gravity field as the distance increases. 

Figures 25 to 27 illustrate this phenomenon with respect to a recent geoid due to W. M. Kaula 
(Kaula-combined [1964] in Table 1). These graphs are a plot of the equatorial longitude force field 



Figure 25— Equatorial longitude gravity forces at 200 miles from the composite geoid of W. M. Kaula (1964). 
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LONGITUDE PERTURBATION (10" 5 g) 



Figure 26— Equatorial longitude gravity forces at 1000 miles from the composite geoid of W. M. Kaula (1964). 





EQUATORIAL, LONGITUDE -GRAVITY ACCELERATION (10‘ 7 g,) 



Figure 28— The earths longitude gravity field on the equator at the synchronous altitude, 

according to geoid (2), (Table 1). 


given by Equation 118 at various altitudes for the longitude coefficients of this geoid. They indicate, 
for example, that orders of longitude gravity higher than the second order can probably be safely 
ignored at almost all longitude locations with respect to their influence on the high altitude 24- hour 
satellite compared to the influence of second-order longitude gravity. On the other hand, they 
probably cannot be safely ignored with respect to their comparative influence on satellites orbiting 
below 6,000 nautical miles. Figure 28 strongly suggests that the recent determination of J 22 from 
24 hour satellite data (References 8 and 11) should be a relatively secure one. 

Figure 29 illustrates the relative imprecision with which the full longitude gravity field of the 
earth at 24 hour altitudes was known prior to the direct measurement by its effect on the Syncom II 
satellite (see References 8, 11 and 15). 

IX. ON THE ORDER OF MAGNITUDE OF THE EARTH GRAVITY PERTURBATION 

DUE TO SMALL WANDERING OF THE SPIN AND AXIS OF FIGURE OF THE EARTH 

In Section VH it was pointed out that there are two consequences from the fact that the earth T s spin 
axis is not quite a principal axis of inertia. In the first place it implies from dynamical considerations 
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Figure 29— Full field earfh longitude gravity forces (perturbations) around the equator 
at 24 hour altitudes (19,300 n.m.) according to three recent geoids. 

that the north pole, the instantaneous spin axis, must wander with a periodic motion and with small 
amplitude both in inertial space and with respect to the body of the earth (Reference 5). This effect 
must be present even in the absence of external torques on the earth, which, of course, can be ex- 
pected to exaggerate it. Such a short-period wandering of the pole (not to be confused with the 
long 26,000-year period of general precession) has in fact been known to astronomers for some 
time (Reference 5, p. 383). The largest component of it appears to have a period of about 440 days 
and an amplitude of about 13 feet with respect to the surface of the earth. The fact that physically, 
the pole is not in the same place on the earth from instant to instant means that there is a small 
periodic latitude variation in the stations used for determining the potential of the earth. Such a 
small variation in latitude cannot mean a significant change in earth potential constants between 
sets reduced from data at different stations at different times or the same station at different times, 
because these constants are by far more sensitive to the radius from the center of the earth than 
any other coordinate. Nevertheless, a 13 foot arc at the earth’s surface is 85 feet at synchronous 
altitudes and within the discriminating power of Doppler range and range rate orbit determination. 
This geometric effect of the consequence of the pole’s short period wandering has not heretofore 


46 



been included in satellite determination programs or in geodetic applications from them. The 
earth constants themselves (as will be shown shortly) cannot be significantly affected by this small 
pole wandering. Putting it another way, physical perturbations, as distinct from the geometric 
effect above, from the slight wobbling of the spin axis will also be shown to be small. Nevertheless, 
the geometric effect appears to be of sufficient magnitude to be determinable from observations of 
satellites. This effect might be discriminated in an orbit determination program which has the 
capability of correcting for the latitude error (on a time basis) in the coordinates of the observation 
stations. Such information on the details of short-period, pole wandering would be of great value 
to geodesists in interpreting the physical nature and behaviour of the interior of the earth. 

We now assess the magnitude of the gravity potential term which arises only when one uses the 
instantaneous spin axis as a coordinate reference and not the true "axis of figure" for the earth (its 
principal axis of inertia). 

Consider in Figure 30 a polar section of the earth containing this true "axis of figure." 

Instantaneously, the axis of figure of the earth X 3 ' is separated by y from the polar axis x 3 . It 
is shown below that the gravity effect of this "tipping" of the reference spin axis is proportional to 
the product of inertia of the earth with respect to the spin axis and the tipped equatorial axis X 2 
(evidently close to zero since y is always very small). 

In Figure 30 note that dm, * , S 3 , <f 2 ' and <f 3 ‘ are all in a plane parallel to the X 3 ' , X 3 plane. 

The product of inertia of the earth with respect to the spin axis, X 3 , and the tipped equatorial axis, 
X 2 , is 


£ 2 £ 3 dm . 

J D 

e 

From the figure it can be seen that 

^2 = ^2 cos y ~ €3 sin ' y 

£3 - €2 sin ^ + £3' cosy (120) 

since the X 2 , X 3 axes are rotated by y in the 
polar section with respect to X 2 , X 3 ' . 

For purposes of simplicity (and with es- 
sential accuracy) the earth is assumed to be a 
biaxial oblate ellipsoid in this argument so that 
every polar section (as in Figure 30) contains 
the second axis of inertia (X a * ) with respect to 
which the moment of inertia of the earth is B . 


X3 



Figure 30— Polar section of the earth containing 
the axis of figure. 


47 


From Equation 120 the product of inertia of the earth with respect to x 2 , x 3 is 


f 2 ^ 3 dm ~ I [^2 cos y ~ £3 ' sin y ) {^2 sin y + €3 cos y ) dm 

j d j D 



+ 


1 


Si S 3 dm 


(cos 2 y - sin 2 y 


} • ( 121 ) 


But since, in this argument, x 2 ‘ is a principal axis, 


I 


f a ' f,' dm = 0; 


and by definition 


and 


J d [(^) 2 + (^) 2 ] dm = c 

U e 

f [fc7) 2 + (^i') 2 ] dm = B • 


Equation 121 therefore reduces to 


f 2 ^3 dm * \ (C - B) s in 2y . (122) 

J D c 

Additionally, it is shown below that the products of inertia with respect to axes x 1? X 2 ( J f j f 2 dm^ 
aildX >’ X 3 (/ D ^1^3 dm ) are both zero. (Not shown in Figure 30 are the and Xj' coordinates of dm 
which are £ x and respectively.) These results hold for this argument because from Equation 120 
and the fact that € x = ^ x9 


£ ! f 2 ~ €2 f i # cos y “ f 3 ' f i' sin y 

^1^3 ” fi' fj sin >' ^ fi' f 3 ' cos y • 
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Since X/, X 2 ', X 3 ' are principal axes of inertia. 


I i ! * " cosy 

j D •'D 


£ 2 dm - cos y f dm - siny I £ 3 £j' dm - 0 


1 


(123) 


and 


J fj^jdm - sinyj 


dm - sin y I dm + cosy 


f Si Si 

•*D 


dm - 0 


(124) 


Other moments of inertia with respect to the X t , X 2 , X 3 axes, important in what follows, can be 
similarly shown to be: 



[(.f,) 2 + (£ 3 ) 2 ] dm = B + (C - B) sin 2 y 

| = B ’ 
e 

[(^l) 2 + (^ 2 ) 2 ] ^ = C'(C-B)sin 2 y 


(125) 

(126) 

(127) 


We now return to the main argument. When the earth was referenced to its principal axes of 
inertia, it was seen that the moment of inertia of the earth about the line to the field point reduced 
simply to Equation 20b or 


However with X,, X 2 , X 3 no longer principal axes, but rotated in the manner of Figure 30 with 
respect to the principal set X,' , X 2 ; , X 3 ' , this line moment of inertia becomes, from Equations 20a 
and 122-127 


r 2 I r = (xj) 2 B + (x 2 ) 2 [B + (C " B) sin 2 y ] + (x 3 ) 2 [c - (C - B) sin 2 y ] - X 2 X 3 (C - B) sin 2y .(128) 

The moment of inertia for the biaxial earth of this simplified argument with respect to its c.m., 

I 0 - i(A + B + C) = ^(2B+C) - \ [2B + c][(x ,) 2 + (x 2 ) 2 +(x 3 ) 2 ]/r 2 . (129) 
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is still given by Equation 20c since it is independent of the coordinate reference. When the argument 
leading to Equation 48 is repeated, but with respect to these altered axes, Equations 128 and 129 in 
(47a) give 





+ ( X 1 X 2) [ 6 ( F 22 )l] + ( X 1 X 3) [ 2 ( F 2 l)o] + ( X 2 X 3) [ 2 ( F 2 l)l] 


(*.)’[ 


2B + C 


3B 

~T 


] * («.)■ . 


2B + C 


3 

2 


{B + (C -B) 


s irr y 


>] 


+ ( x 3 ) 2 [ 2B 2 C " “ y {c - ( c - B) sin 2 y \ [x 2 x 3 (C - B) sin 2y] . (130) 

It is to be understood that the longitude reference (the X x axis) has not changed in the particular 
pole tipping of Figure 30. Since the earth is almost a biaxial ellipsoid, one can always choose the 
longitude reference axis arbitrarily or describe the tipping in this particular fashion. The physical 
effects which follow should be the same as for an arbitrary pole migration except for a longitude phase 
shift which would depend on the longitude as well as the latitude variation of the pole. Such a 
variation could be described with respect to the Greenwich meridian through the axis of figure. 

The gravity constants ( F 20 ) 0 > ( F 22 )o > ( F 2 2)17 ( f 2 i)o> ( F 2 1 ) 1 in Equation 130 thus refer to axes 
Xj, X 2 , x 3 as in the argument of Equation 48. But now, since these axes cannot all be principal 
axes, we have chosen X t and x 2 arbitrarily to simplify the description of the effect. 

Equating coefficients of linearly independent terms on the right and left sides of Equation 130 
as in Equations 49 - 54 gives the following: 


( F 2 0)0 . 1 

— 2— + 3( F 2 2 )o = 2 (C-®) 

" ( F 2 0)0 . C-B 3 1 / , \ 

2 3 ( F 22 ) 0 “ — 2~ - 2 Sln2 7 (C - B) - 2 ( c _ B > V 1 ' 3 Slr|2 y ) 


( F 2 0 ) 0 = B - c + y [sin* r (C-B)] - (B-C) (l - f sin 2 yj 

6( F 22), = o 

^( F 2l)o _ 0 

2 ( F 2 i ) j = T (C-B) sin 2y . 


(131) 

(132) 

(133) 

(134) 

(135) 

(136) 
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- From Equations 131 - 136 the gravity constants in terms of the moments of inertia with respect 
to the axis of figure and the tip angle y are 


- (B - C) (l - y sin 2 y) 

(137) 

, - \ (C -B) sin 2 y 

(138) 

(*ai)i = 0 

(139) 

K) 

O 

II 

O 

(140) 

= f (C -B) sin 2y 

(141) 


Note that when y = 0, the results of Equations 137 - 141 are the same as Equations 49 - 54 
with A = b. In other words, the effect of a small tipping of the axis of figure with respect to the 
spin axis is to change j 20 by the order of magnitude of sin 2 y, to change j 22 by the order of magni- 
tude of j 20 sin 2 y , and to introduce an additional J 21 longitude effect. All of these gravity effects 
are very small and in most applications should be negligible, as will be seen below. 

When Equations 43 and 56 are rewritten with respect to the non-principal axis set X„ X 2 , X 3 , 


F [F 3F 1 

v e - — — + — |~ 2 - (3 sin 2 L - l) + — 3"~ sin 2Lcos ($~<£ 2 ,) + 3F 22 cos 2 Lcos 2($-<f> 22 ) 1 , (142) 

and 


GM j J 

V e = — — 'll j- (a/r) 2 (3 sin 2 L - l)~ -y (a/r) 2 J 21 sin 2Lcos <j> 2 j) 

- 3J 22 (a/r) 2 cos 2 Lcos 2 (<£ _ <£ 22 ) ^ • (143) 

For the J 21 , f 21 coefficients, the identification of Equations 142 and 143 implies 

3F 21 _ -3(i/r) 2 J 21 CM e 

2r 


or 


-F 


J21 


21 




GM e (a) 


(144) 
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Similarly, the identification of Equations 142 and 143 gives 


J 20 


-F 


20 


GM e (a) 2 


(145) 


Equation 140 in (46a) implies <p 21 - n/2 or 2m / 2 radians, which implies F 21 - 2G(f 21 ) 1 ^3 . This 
result and that of Equation 141 give 


or from Equation 144 


F 21 » f [|(C-B)]sin2r 


(C - B) sin 2y 


(146) 


Since F 20 = G(F 20 ) 0 , (see Equation 46a), from Equations 145 andl37, j 20 is given approximately by 


C -B 
M e (a) 2 


(147) 


Combining Equations 146 and 147 gives 


'21 


I 2 0 


sin 2y •_ 
2 


* 20 


(148) 


since y«l. But 


y ~ (3960x 5280)' 0.62 x 10 6 , 

(Reference 6, p 383) and J 20 = 1082.5 X 10~ 6 . Therefore, from Equation 148 

|j 21 | = 0^00067 x 10“ 6 . (149) 

Earth longitude coefficients |j nm | (from Equation 99 for example) determined by ignoring J 21 
altogether, all appear to be greater than 0.01 x 10~ 6 for n < 6 (Reference 12). It is clear 
then that for orbit determinations which need inclusion of earth tesserals up to about the sixth 
order, the physical effects of the small periodic wandering of the north pole are entirely negligible. 
In fact, if J 21 is arbitrarily set to zero, it in effect forces the latitude reference to be the true axis 
of figure at any instant in an orbit determination which utilizes an earth potential such as Equa- 
tion 99. In the torque-free nutations of the oblate earth the spin axis is always considerably closer 
to the "invariable” angular momentum axis (the true inertial reference for motion with respect to 
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the earth) than the axis of figure. Therefore, the time based correction for station latitude and 
longitude bias from observations of satellite motions will, assuming torque-free nutations only, 
give almost directly the movement of the axis of figure in inertial space. 
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Appendix A 

Notes on Notation 


It would appear that almost every author of geodetic investigations uses a different notation or 
definition for the constants in the gravity potential of the earth, for example, Reference A-l. 

W. M. Kaula (Reference A-2) has reviewed these and has proposed a limited number of standard 
forms useful in the applications and easily converted from one to the other. They are all based on 
the form of the potential as it is written in Equation 99. As can readily be appreciated, it is 
important that when using a potential in this form, the earth radius (a, for example) be clearly 
defined and stated when constants are reported. 

In one form the longitude constants (m / 0) are 

K = -j (A-l) 

as written in Equation 99. By definition the J nm (m ^ 0) in Equation 99 are all negative. With this 
definition the K nm are interpreted physically as in Figures 5 to 15 for both the K nro and J nni form for 
the potential. There is no logical reason why one of these forms should not be dropped. Historically 
(Reference A-3, for example) the J nn! form of Equation 99 has always been used for the zonal 
(m = 0) potential. It has only been fairly recently, to the author's knowledge, that the K nm form 
(with K nm > 0 by definition) has been adopted for just the longitude part of the series. The author 
feels that the form of the J nin series valid for all n and m in Equation 99 is a reasonable unifying 
compromise of the historically different forms (j n and K nm ) for treating zonal and longitude gravity. 
Again it is noted that the J nm 's in the form of Equation 99 must always be negative so that the k nm 
in the J nm form refers to the same longitudinal axis of symmetry as in the K nm form for the nm 
harmonic. For example, if J 22 in the potential series 99 were positive, this would mean that X 22 
would locate the minor axis of the elliptical equator instead of the major axis as it does when J 22 
is negative. In the K nm form of the potential it has been conventional to define K nm always positive 
so that k 22 locates the major equatorial axis. The simple assignment, J nm = -K nn , preserves the 
same physical interpretation of the k nm between the two forms (Figures 5 to 15). 

The third standard form proposed by Kaula is a simple expansion of the K nm form for 
longitude-dependent gravity. When the longitude -dependent part of the general longitude gravity term 
is written as 

K nra cosm(\-X nn ) ° (k nm cosmXj cos m\ + (k nm sinm\J sinm\ , 

the two longitude constants for the nm harmonic can be defined as 

C = K cosm\ „ 

n rn nm mn 
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and 


S n„ = K nm sinnA nm (A-2) 

With these new constants the harmonic potential series of Equation 99 is (Reference A-2). 


n 1 00 n 


(a/r ) n P n m (sin L) (C nm cos mX. + S nm sinm\) 


(A-3) 


Alternately the K , \ form is 

^ nm ' n m 


. . , 00 n 

% |. \ 


F j 1 ( ^/r ) n P n m ( s in L) K n m cos m(x. — m ) 


n = 1 m = 0 


(A-4) 


and the j „ . k m form is 

n m ' n m 


e r 


co n 

l-££ ( i/r ) np n m (sioL) Jnn, cosm (^“\,rn) 


n - 1 m - 0 


(A-5) 


The C nm , s nm form in Equation A-2 is often a convenient calculating form as the effects of the two 
longitude constants are linearly separated in it. 

From Equations A-l and A-2 the connections between the K , \ , j , C , and S nm constants 

* nm ’ nm" ‘'nm* nm* nm 

are summarized as 


Jn m 


-K 


n m 


J = -(c 2 +s 2 ) 1/2 , 

J n m \nm nm/ 

Km = -S’ tan_1 (W C n m ) • 


and 


• (c 2 + s 2 ) 

\ nm nm/ 


1/2 


(A-6) 


In addition to Equations A-3 to A-5 Kaula (Reference A-2) proposes a form of Equation A-3 which 
defines the potential as 


V 


e 



+ 


00 n 

( i ' / O np n"’ (s inL )[C nm cosnA + S nm 

n = 1 m “ 0 



(A-7) 
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where 


( ^nm » ^nm ) 


(n + m) ! 

(n - m) ! ( 2n + 1 ) (2 - 8,, 0 ) 


1/2 


(c , s ) 

\ n m n m / 


with 


5 m ° - 1 for m H 0 (zonal terms) 


and 


& m ° - 0 for m t 0 (longitude terms) . (A-8) 

It can be shown (Reference A-4) that C nro and S nm are coefficients of harmonics which have a mean 
square amplitude of 1 over the geoid for all values of m and n. Thus, comparison of the order of 
magnitude of the physical effects of the harmonics at any altitude is possible using (C nm , s nm ) values 
modified by (a/r) n to account for the decrease in harmonic amplitude with distance from the c.m. 
of the earth. 


REFERENCES 

A-l. Townsend, G. E., "Handbook of Orbital night, Vol. 1, Part 1," NASA SP-33, 1963. 

A-2. Kaula, W. M., "A Review of Geodetic Parameters," Goddard Space Flight Center Document 
X-640-63-55, 1963. 

A-3. Kozai, Y., "Numerical Results from Orbits," Smithsonian Inst. Astrophys. Obs. Spec. Rept. 
101:1-12, July 31, 1962. 

A-4. Kaula, W. M., "Tesseral Harmonics of the Gravitational Field and Geodetic Datum Shifts 

Derived from Camera Observations of Satellites," J. Geophys. Res. 68(2):473-484, January 1963. 


57 


Appendix B 

List of Symbols 




Qi 

ds i 

r , L, \(or 4>) 

V 

R, H, $ 
D m 

n 

pm 

n 

P n°° r P n 

F n m 

( 0r< ^„ m ) 


G 

M 


x. or x.' 

i i 


A, B, C 

J n m 


differential operator signifying d 2 /** 2 +d 2 /dy 2 +d 2 /dz 2 in x, y, z rectangular coordinates, 
for example 

metric coefficient of the coordinate q. (dimensionless) 
general orthogonal coordinate 

small line element in the direction of increasing q. only 

spherical coordinates of radius, latitude and longitude with respect to the mass distribution 
gravitational potential (units of energy/mass) 

separated variations of the potential; functions of radius, latitude and longitude, respectively 
intermediate potential constants (units of energy/mass) 
associated Legendre function of order n, power m 
Legendre polynomial of order n 

intermediate potential constants (units of energy/mass) 

longitude of the principal plane of symmetry corresponding to the nm harmonic of the 
potential of the mass distribution 

universal gravitational constant (6.673 X 10" 8 dyne-cm 2 /gm 2 ) 

total mass of the distribution whose potential is desired (M e = M = Mass of the earth, in 
the applications in this report). 

moment of inertia of the distribution about the origin 

moment of inertia of the distribution about the line vector 7 to the test point in the field 

rectangular coordinates in inertial space aligned to the axes of principal moment of 
inertia of the mass distribution 

principal moments of inertia of the mass distribution 

potential constant for the nm harmonic of the potential of the mass distribution (units of 
energy/ mass) 
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a, R 0 average of major and minor radii of the elliptical equator of the model earth 

f polar flattening or oblateness coefficient of the model earth ellipsoid (dimensionless) 
e eccentricity of the elliptical equator of the model earth 
a, b, c halt axes of the model earth ellipsoid (c is the polar radius) 
r radius to the surface of the model earth ellipsoid 
^ rotation rate of the earth (units of angular velocity) 

c 0 constant amount of the potential at the surface of the model earth ellipsoid 
a centrifugal potential constant of the earth 

\ geographical longitude measured east from the Greenwich meridian 

A r deviation in radius from the surface of the average earth sphere of radius a to an 
equipotential surface 

K 0 , K 0 intermediate potential constants 

D domain of the mass distribution 

p distance from the field point to a point mass of a mass distribution 

€ the distance from the origin of coordinates to a mass point in the domain of the mass dis- 
tribution whose gravity potential is sought. 
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